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We report the development of a new ray-tracing simulation tool having the potential of the full 
characterization of a radio hnk through the accurate study of the propagation path of the signal from 
the transmitting to the receiving antennas across a perturbed atmosphere. The ray-tracing equations 
are solved, with controlled accuracy, in three dimensions (3D) and the propagation characteristics 
are obtained using various refractive index models. The launching of the rays, the atmospheric 
medium and its disturbances are characterized in 3D. The novelty in the approach stems from the 
use of special numerical techniques dealing with so called stiff differential equations without which 
no solution of the ray-tracing equations is possible. Starting with a given launching angle, the 
solution consists of the ray trajectory, the propagation time information at each point of the path, 
the beam spreading, the transmitted (resp. received) power taking account of the radiation pattern 
and orientation of the antennas and finally, the polarization state of the beam. Some previously 
known results are presented for comparative purposes and new results are presented as well as some 
of the capabilities of the software. 



I. INTRODUCTION 

Multipath propagation is believed to be the major 
cause of data transmission impairments in terrestrial line 
of sight microwave radio systems. Efficient antenna de- 
sign requires the understanding of the propagation of in- 
dividual rays across the channel and gauging the refrac- 
tive index of the various atmospheric disturbances any 
given ray encounters during its propagation. Adopting 
a refractive index model for a given disturbance arising 
from spatial fluctuations in humidity, pressure or tem- 
perature (these fluctuations might be temporal as well, 
but we shall consider, for the time being, that the prop- 
agation time occurs on a time scale much smaller than 
the one associated with these fluctuations), we establish 
the ray propagation equations and solve them with sev- 
eral numerical techniques having a first, fourth and sixth 
order accuracy. The ray tracing equations are initially 
solved in two dimensions bypassing the effects of small 
and non-linear terms as explained in section 2. Later on, 
we switch to 3D in order to assess the effects the small 
and non-linear terms have on ray propagation. Several 
facts emerge from this approach: 

• The small non-linear terms lead to a breakdown of 
standard integration techniques. The ray equations 
which constitute a system of 6 ordinary coupled 
non-linear differential equations become stiff. This 
means the integration step becomes so small (be- 
cause of the presence of terms that differ by several 
orders of magnitude) making the integration pro- 



cess so slow that any progress in seeking a solution 
of the system is virtually stopped. 

• The relation between the launching and arrival an- 
gles for a given disturbance are profoundly altered. 
What was previously believed to be a "good" or 
"bad" launching angle might have gotten its true 
attributes from reasons different from what is cur- 
rently known. 

• A very high sensitivity is observed around certain 
launching angles: a very small uncertainty in the 
launching angle can induce the ray to take a path 
radically different from what is normally expected. 

This report is organized in the following way: In sec- 
tion 2, we establish the ray-tracing equations (RTE). In 
section 3 we describe some of the problems encountered 
during the solution of the RTE, namely those related to 
stiffness and present the algorithms to cure them (Ap- 
pendix A contains a description and an example of a 
stiff system). In section 4 we compare our approach to 
previous ones and present some illustrative new cases in 
section 5. This section also describes the potential ap- 
plications of the software and its capabilities. Section 6 
discusses some possibilities for future developments. Ap- 
pendix B shows how to avoid stiff differential equations in 
two dimensions and turn the RTE into a set of recursion 
relations. 



II. RAY TRACING EQUATIONS 
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In terrestrial microwave radio systems, the range of 
frequencies used and in comparison the range of length 
scales present in the channel allow us to use a geometric 
(or ray) approach to electromagnetic propagation. The 
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fundamental equation of geometrical optics is the Eikonal 
equation : 



(gradS*) 



(1) 



where n is the local refractive index and S is the local 
phase of the ray. Taking the gradient of both sides of the 
Eikonal equation gives the second order vector propaga- 
tion equation: 



ds 



gradn 



(2) 



where R is the ray position and ds is a differential 
displacement along the ray path, i.e. ds = ||dR||, the 
norm of the vector dR. 

This can be rewritten as a system of two first order 
equations: 



dR 

ds 
d{nT) 
ds 



T 

gradn 



(3) 



where T is a unit vector tangent to the ray path (The 
geometry is depicted in Fig.l). The advantage of solving 
a first order system rather than a single second order 
system is threefold: 



• Stability problems are easier to handle. 

• Validity of the solution is easy to monitor since one 
has to have for all times ||T|| = 1 providing a sim- 
ple means to check the quality of the integration 
procedure. 

• Accuracy of the solution is controlled within certain 
tolerance limits depending on the selected integra- 
tion step. 

This is discussed in detail in section 5. The refractive 
index function of the atmosphere is written as: 



where k is the refractive index gradient with height 
h. The atan() term above is due to a disturbance lo- 
cated at a height ho having an extent Ah and a refrac- 
tive strength An. For a normal atmosphere (An = in 
the Webster model) both models are linear in h (after 
expanding the exponential to first order). Nevertheless, 
their dependence solely on height does not account for 
the 3D nature of the atmosphere and its disturbances. 
Some models like the recent one introduced by Costa ||] 
mimics a 3D atmospheric disturbance by multiplying the 
refractive index along the vertical with a Gaussian func- 
tion along the horizontal perpendicular to the ray path 
plane. Going beyond these approaches, we introduce a 
full 3D profile: 



N = p^{x)py{y)ph(h) + kh + No 



(5) 



where Px , Py and ph are the index profiles of the dis- 
turbance along the three directions in space x, y and h. 
A^o is an average normal atmosphere index and k is the 
index gradient along the height. A profile function p{X), 
along direction X is typically taken as: 

p{X) = {Anx/2)[tanh{{X - Xi)/AXi) 

- tanh{{X - X2)/AX2)] (6) 

where Xi (resp. X2) is the point where the hump 
starts growing (resp. decaying) and AAi (resp. AA2) is 
a typical length scale for the growth (resp. decay). An^; 
is the refractive strength of the disturbance. This model, 
though realistically representing a localized anisotropic 
disturbance in the atmosphere is based on a separable 
model of the refractive index function. 

While our methodology can handle any arbitrary 3D 
model of the refractive index, any of these refractive mod- 
els have to be modified in order to take account of the 
curvature of the Earth by the inclusion of a term |^ equal 
to lO^h/Re where Re is the radius of the Earth. 



III. STIFF DIFFERENTIAL EQUATIONS 
ALGORITHMS 



n = 1 + IQ-^N 



(4) 



where N depends on the frequency used, humidity 
conditions and height above the Earth ground. Several 
models exist for the range of frequencies and heights we 
are dealing with and are generally expressed in N units. 
The following two models are of interest; the first for a 
normal atmosphere and the second for a disturbed one: 



• Exponential model: N = 315 exp(— 0.136/i), with 
h (height) in kms. 



• Webster 

^ atan(12.63 

7r ^ 



model: N 

ih~ho)_\ 
Ah I 



300. 



kh 



Using [4] , the ray-tracing system [3] is rewritten as: 



dR 

ds 
dT 

ds 



= [gradA^-T(gradAr.T)]/(A^+ 10^) (7) 



Two important features appear in the RHS of the sec- 
ond equation in the system: 

• The non- linear term in T. 

• The wide range of orders of magnitudes in the de- 
nominator. 
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These terms can be eliminated with the following pro- 
cedure: Replace equation [7-b] by another equation defin- 
ing the curvature of the ray path r: 

f=U/p (8) 

where U is the normal to the trajectory. U is perpen- 
dicular to T and normalized: ||U|| — 1. The unknown p 
can be determined by taking the scalar product of both 
sides of [7-b] with U and using [8]; one gets: 

l/p==U.gradiV/(7V + 10^) (9) 
Substituting [9] in [8] gives the following system: 



— = U(U.gradiV)/(A^+ 10^) (10) 
as 

In general, this system is not closed because it involves 
U besides R and T. In two dimensions, one can close 
the system by invoking [|| the orthogonality of U and T 
through: 

U = xxT (11) 

where x is the unit vector along the x direction. With 
relation [11], system [10] is now closed and can be in- 
tegrated by any standard explicit integration method 
(Predictor-corrector, Euler, Runge-Kutta, Richardson 
etc.). This will be illustrated in section 4. In gen- 
eral, iV is a function of the position vector R; when it 
depends only on the height, it is possible to further sim- 
plify the system and reduce it to a single scalar equation. 
In the case N depends only on height, gradA^ is along 
the vertical and if ij) is the angle T makes with the lo- 
cal horizontal, U being perpendicular to T will make the 
same angle with the vertical, [9] yields: 

- = \{dN/dh)cos'iP\/{N + 10^) (12) 
P 

Livingston [Q has derived an equation similar to [12]: 

- = -{l/n){dn/dh)cosil; (13) 
P 

Equation [13] is equivalent to [12] when the right 
sign is used. We have integrated system [10] in two 
dimensions and recovered typical results found in the 
literature, avoiding the difficulty arising from [7-b]. In 
the three dimensional case, one has to deal directly 
with system [7] with all terms retained, for, in general, 
the T vector does no longer have to be confined to the 



transmitter (TX) receiver (RX) plane. In this case, all 
standard explicit integration schemes break down. In 
other words, the norm of the vector T tangent to the 
ray path is no longer conserved. In order to fulfill the 
condition ||R|| = 1, one has to take an integration step 
so small that the integration process is virtually stopped. 
This is called stiffness and an illustrative example is 
given in Appendix A. 

Stiffness can be cured with the so called implicit in- 
tegration schemes. In contrast to explicit integration 
schemes where a current system value depends only on 
the previous ones, implicit schemes couple present and 
past values of the system altogether. A price to pay is an 
increase in CPU time but the rewards are stability, accu- 
racy and large integration steps. We have implemented 
two implicit schemes: 

• Generalized Runge-Kutta (GRK) method of fourth 
order [^. 

• Roscnbrock (ROW) method of sixth order [^ . 

In the first scheme, given a system of first order ordi- 
nary differential equations (ODE): 

dy/ds = f (y) (14) 

one builds the vectors from the system values at step 
n-1: 

ki = erf (y„-i + ^ aijkj) with: i, j=l...m 

(15) 

and evaluates the next value n of the system with: 

Yn = Yn-l + ^ hk, (16) 

a is the integration step and the and bi are coef- 
ficients depending on the scheme m of the integration 
order. In the Roscnbrock case, one adds to [15] the 
term cr(^)^(iyki, where the d^^s are order dependent 

coefficients and (^) is the Jacobian of the system. The 
above equations are implicit since the unknown vectors 
ki needed for integration step n appear on both sides 
of [15]. In the GRK method, only the vector function f 
is needed whereas in the ROW case both f and its first 
order derivative (Jacobian) are needed. 

Both methods have been proven to perform very well 
up to stiffness parameters (ratio of the highest to the 
smallest eigenvalue of the Jacobian) as high as 10*". Inci- 
dentally, our stiffness parameter has been observed (while 
testing ROW algorithms) to be generally around 10*. We 
have used GRK of order 4 and ROW of order 6 because 
they have been extensively tested for a wide range of 
systems and are thoroughly documented. 
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IV. VALIDATION OF THE APPROACH AND 
COMPARISONS WITH PREVIOUS 
TREATMENTS 

In order to validate our technique, we started with a 
comparison against analytically known solutions. Three 
models were tested, the axial gradient refractive index 
case, the sine- wave optical paths and the classical Luneb- 
urg lens (see, for instance, reference 7). In all three cases 
our results compared very accurately with the analyti- 
cal ones. Then we went ahead and proceeded to solve in 
detail a case well documented in the literature and inves- 
tigated by Webster |^ for various launching angles. This 
model is two dimensional (2D) and extensively referred 
to in the literature. We use the 2D version of the system 
of equations [10] which is non-linear (N is a non-linear 
function of R and a power of U appears in [10-b]). 

The integration, started by taking values of R and 
T as the initial location and launching vectors, is done 
with a first- order Euler and fourth order Runge-Kutta 
methods. The TX-RX configuration and propagation 
conditions are the same as those given in Table 1 of 
Webster's @] paper. In Fig. 2 we show the various 
ray paths between the TX and the RX for a series of 
launching angles (taken with respect to the horizontal) 
varying from -0.25 up to 0.5 degrees. The different 
launching angles, we use, are respectively, in degrees: 
-0.25, -0.20, -0.15, -0.10, -0.05, 0.0, 0.10, 0.20, 0.30, 0.40, 
0.50. The refractive index profile used in the study is 
displayed in Fig. 3. 

While Fig. 2 is based on a first order (Euler) integra- 
tion method, some changes might occur if we rather use 
a fourth order Runge-Kutta method. In fact, the ray 
paths based on either scheme show no appreciable dif- 
ferences and compare well with the results found earlier 
by Webster in the same conditions. However, some dis- 
crepancies appear for positive launching angles and are 
probably due to the different levels of numerical accu- 
racy between our treatment and Webster's. Let us recall 
that in our case the numerical accuracy is monitored by 
checking the conservation of the norm of T. In these 
simulations, it is conserved with an error smaller than 
10^^. In order to compare our results to Webster's di- 
rectly, we derive, in the same fashion, recursion equations 
for the ray radial distance R (taken from the center of 
the Earth) and the angle 4' that T makes with the local 
horizontal. Referring to Appendix B and Fig. 4, we can 
write the following relations: 

i?2 = Ri + ds sin[ipi) (17) 

= ^IJi+ds — ^P^-sm i( — ) (18) 
Ri pi 

where the radius of curvature pi is given by [12] with 
■0 = ^Ai and dN/dh is taken at the height R — R^ (R^ 
is the Earth radius). For a given step ds, one starts 



the set of iterations [17] and [18] with the launching 
radial distance Ri and angle ipi. Using the same 
initial values as before we retrieve almost the same 
ray trajectories obtained in Fig. 2. The validity of our 
results is monitored by the constancy of the modulus 
of T versus 1. Additionally, we compared our results 
(Euler and Runge-Kutta) to a very high accuracy inte- 
gration technique based on the Butcher's [^ algorithm 
(seven-stage sixth-order Runge-Kutta scheme). The 
sixth order results are virtually identical to the fourth 
order's and Fig. 5 depicts the ray trajectory obtained 
with the different levels of accuracy under the same 
atmospheric and launching conditions. Incidentally, 
the difference between fourth and sixth order trajecto- 
ries in Fig. 5 are on the order of a fraction of a millimeter. 

In spite of the above agreement, which is basically rela- 
tive, one still has to gauge independently the accuracy of 
the results for a selected order and integration step. This 
is done with the following method: Pick an order p and 
an integration step a; integrate once with a and twice 
with cr/2 in order to reach the same point; define a step 
ratio K from the difference A between the two results: 



K= "+^2^/(2^-^ l)(A/e) (19) 

and monitor the value of k for a given tolerance, during 
integration. Ideally, we should have k < 2. In Fig. 6, we 
display k versus the integration step number for the first 
order (Euler, p=l) case as well as the Runge-Kutta 4-th 
order (p=4) and Butcher 6-th order (p=6) for a tolerance 
of 1 millimeter. We use exactly the same condition as 
previously and a launching angle of 0.2 degrees. The 
figure shows clearly the superiority of 4-th and sixth order 
methods for the selected step when such a high accuracy 
is desired. 



V. ILLUSTRATIVE RESULTS AND 
CAPABILITIES OF THE METHODOLOGY 

We move on to the description of the 3D propagation 
case and show, with a simple example, how we evalu- 
ate the power from the antenna radiation pattern, the 
beam spreading and the state of polarization. We select 
a coordinate system such that the TX is somewhere on 
the z-axis whereas the y-axis is along the TX-RX line. 
The vertical plane is defined by the z axis and the TX- 
RX line. The beam spreading is evaluated by launching 
simultaneously several beams in the vertical and horizon- 
tal planes with angles differing by a small amount from 
those characterizing the main beam. The logarithm of 
the ratio of the surfaces swept by the different beams at 
the receiver location gives an estimate of the spreading 
loss. In order to account for the TX-RX antenna ra- 
diation pattern, we simply recall that the electric field 
radiated by a parabolic circular aperture antenna at a 
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point defined by its distance r from the main lobe origin 
and making an angle 9 with the lobe axis is given by: 



E(r,e) ^ jPEoa[exp{-jPr)/r] Ji{(3a sind)/Psind (20) 

where a is the aperture radius, Eq is a reference field, 
f3 = 2n/\ with A the wavelength used, Ji is the Bessel 
function of the first kind and j = y/—T. The antenna pat- 
tern is obtained after normalizing the value of \E(r, 9)\: 



f{0) = (2//3a)|Ji(/3a sin9)/sin9\ 



(21) 



Alluding to our choice of axes, if the main lobe is point- 
ing in a direction defined by the angles /3, 7 (in the verti- 
cal and horizontal plane respectively) and we have a ray 
along /3', 7', the angle the ray makes with the main lobe 
axis is: 



9 = cos ^{cosP sin^ cos(3' sin'f' 

+ cosP cosj cos/S' cosj' + sinP sinj3'){22) 

The power (in dB) is given by 20 logiof{9). The polar- 
ization state of a ray rotates, during propagation, by an 
angle calculated with the help of the following formula: 



0(^,5) 



ds/r 



(23) 



where A and B represent the two end points of the 
ray trajectory; r, the local torsion of the ray is different 
from zero when the trajectory is not confined to a plane. 
Using the Frenet-Serret jl] formula: 



dB/ds 



-rU 



(24) 



Taking the dot product with U on both sides of equa- 
tion [24] and replacing the value of t in [23], one gets: 



(j){A,B) = dsV(dB.U) 

J A 



(25) 



In order to evaluate the polarization rotation of the 
ray propagating from Ato B with [25], a finite difference 
approximation B„ — B„_i is used for the differential dB, 
where the subscripts refer to the integration step. The 
final discrete formula for the polarization angle reads: 



n=N 



(26) 



n=l 



where N is the number of integration steps between A 
and B. For illustration, we treat two 3D examples. In 
the first case, we take a refractive index model consisting 
of a refractive layer of finite length along the TX-RX 



line. The linear extent of the layer is taken respectively 
as 5, 10 15, 20 and 25 kms. Fig. 7 shows the dramatic 
effect the extent has on the ray path. Incidentally, the 
refractive index model along the height is taken as the 
same Webster model as before and the ray launching is 
made in the vertical plane. In the second case, we take a 
refractive index model given by a Webster profile along 
z and a profile Py{y) given by [5]. Moreover we take 
an arbitrary 3D launching direction. The resulting 3D 
ray trajectory for the selected parameters listed in the 
corresponding caption is displayed in Fig. 8. 



VI. CONCLUSIONS AND FUTURE 
DEVELOPMENTS 

We intend to use this technique to study the dynam- 
ics of microwave radio signals controlled by unstable 
atmospheric layers. The instabilities cause short error 
bursts lasting from many tens of micro-seconds to a 
few milliseconds 10 . Since, the error bursts have 
detrimental impact on communication networks [jll] , the 
future digital radio systems should be made immune to 
radio propagation degradations causing them. In order 
to develop defense strategies against the error bursts 
caused by atmospheric propagation instabilities, the 
physical characteristics of the instabilities have to be 
well understood. This 3D ray- tracing technique will 
be used to study the effects of dynamically changing 
atmospheric layers of limited size on microwave radio 
signals received simultaneously by a few parabolic 
antennas [l2| . A propagation model simulating the 
recorded dynamics of received radio signals |lO[ will, 
not only, help understanding the physical causes of the 
error bursts, but it will also be used in the computer 
optimization of antenna designs capable of minimizing 
the frequency of occurrence of the propagation caused 
error bursts. Highly accurate numerical techniques are 
required since small fluctuations of the atmospheric 
conditions are believed to be responsible for the flat 
phase fluctuations impairing the digital demodulation of 
the received microwave radio signals. 



APPENDIX A 



Let us consider the following first order system con- 
sisting of a pair of linear ordinary differential equations: 

dyi/dx = X+yi + X^y2 (27) 
dy2/dx = X-yi + X+y2 (28) 

where A+ = (Ai + A2)/2, A_ = (Ai- A2)/2 and a; > 0. 
The solution of the system is: 
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yi = Ciexp{Xix) + C2exp{X2x) (29) 
2/2 = Ciexp{Xix) - C2exp{\2x) (30) 

where Ci and C2 are constants determined by the ini- 
tial condition at x=0. In order to conform to our no- 
tation of Section 3, we define a column vector y whose 
components are y\ , y2 and write the system as: 

dy/dt = f (y) (31) 
The eigenvalues of the Jacobian of the system: 



R2 = Ri + ds sin{il}\) 



ipi + ds- 



(34) 

sin-\—) (35) 
Pi 



R2 = Ri + ds sin{tpi) 



(36) 



where ds = ||R2 — Ri|| . The radial distance i?i (resp. 
R2) is taken from the center of the Earth. The angle 
56 between the two radial directions may be found by 
inspection: 



di_ 
dy 



A+A_ 
A_A+ 



are solutions of: 



rfei|AI-(^)|=0 



(32) 



(33) 



where I is the (2x2) unit matrix; they are nothing else 
than Ai and A2. If one picks Ai = — 1,A2 = —1000. and 
chooses an explicit integration method, one finds the 
integration step should be smaller than 2/IA2I , which 
is 0.002. This is the origin of stiffness: even though the 
term exp(-1000 x) contributes almost nothing to the 
solution for a; > 0, its presence alone, virtually stops the 
integration process. 



APPENDIX B 



Risin{56) = ds cos(V'i) 
which can be approximated by: 



56 = ds cos{'tpi)/Ri 



(37) 



(38) 



In order to find the relation between the angles 
ipi and ip2, we use the relation defining the derivative 
of T, dT/ds = U//9 in a discrete form: 



T2 - Ti = ds Ui/pi 



(39) 



Taking the scalar product with Ui on both sides of 
above, one gets: 



T2.U1 = ds/pi 



(40) 



The inspection of Fig. 4 provides the angle between 
T2 and Ui: 



The geometry of propagation is shown in Fig. 4. At 
any point along the ray trajectory the tangent vector T 
makes the angle with the local horizontal. When the 
ray propagates between two nearby locations, one may 
write: 



(T2,Ui)=V'2-V'i-'5^ + 7r/2 (41) 
Using the above result gives the relation sought: 

1P2 = tpi+ 56 - sin-\ds/pi) (42) 
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Figure Captions 



Fig. 1: Geometry of the system showing the coordinate 
system, the antennas in the vertical yz plane, a 
typical ray path and the local Frenet-Serret system 
(T, U, B) attached to a point along the path. 

Fig. 2: Euler first order 2D results. The rays are launched 
in the vertical plane and the angle they make with 
respect to the horizontal xy plane is respectively: - 
0.25, - 0.20, -0.15, -0.10, -0.05, 0.0, 0.10, 0.20, 0.30, 
0.40, 0.50 degrees. Equations [10] are used along 
with model [4-b] for a perturbed atmosphere N — 

300. + kh+ ^ atan(12.63^^^^) with the same 
parameters as those given in Table 1 of reference 
2: k = -39, An = -20 (both in N units), /io=175 
meters, Ah —100 meters, the transmitter height is 
125 meters and the TX-RX separation is 60 kms. 

Fig. 3: Webster Q model refractive index function (in 
N units) along the vertical showing an anomaly 
at a height of 175 meters and whose width is 
equal to 100 meters. The curvature of the Earth 
term lO^h/Re is present. The negative gradient 
of the layer refractive index is responsible for the 
multipath effects observed. 

Fig. 4: Geometry of the ray trajectory used for establish- 
ing the recursion equations. The local tangent T 
vectors are shown making the angle ip with the 
local horizontal perpendicular to the ray vectors 
R drawn from the center of the Earth O. Two 
neighboring points along the ray paths are shown. 

Fig. 5: Comparative study of the ray trajectories obtained 
from the recursion relations [17] and [18] (upper- 
most long dashed curve) and 1st order Euler (full 
line curve) on one hand, and the Runge-Kutta 
(4-th order) and Butcher (6-th order) on the other 
(short dashed curve). The launching angle in all 



cases is 0.2 degrees in the vertical plane and the 
model considered for the refracting layer is the 
same as Figure 2. The fourth and sixth order 
results are virtually identical. 



Fig. 6: Comparative study of the behavior of the step 
ratio versus step number for the Euler (1st order), 
the Runge-Kutta (4-th order) and the Butcher 
(6-th order) methods when the step is fixed to its 
starting value. Ideally, this ratio should always 
be about 2. In the first order case, the bound is 
violated very rapidly (upper curve), whereas it is 
respected until almost the end of the trajectory 
in the 4-th (long dashed curve) and 6-th order 
(short dashed curve) cases. The tolerance is 1 mm 
and the step used is one hundredth the TX-RX 
distance. 



Fig. 7: GRK (Implicit, 4-th order, 3D) results for the ray 
trajectories when the extent of the layer is a vari- 
able. Starting with a launching angle of 0.2 degrees 
in the vertical plane, the layer spans, initially, the 
entire hop of 60 kms (lowest curve). Moving up- 
ward from the next lower curve, the layer extent 
(along the TX-RX line) is from 5 to 25 kms, then 
5 to 20 kms, 5 to 15 kms and finally 5 to 10 kms. 
In all cases, the refracting layer model is the same 
as in Figure 2. 

Fig. 8: GRK (Imphcit, 4-th order, 3D) results for the 
ray trajectory when the refractive index of the 
layer varies along two spatial directions (y and z) 
and round Earth profile considered. The normal 
atmosphere parameters are No=SOO N units and 
the gradient k=-39 N units/km. The 3D refractive 
index layer is described with a profile along y given 
by [tanh{{y - yi)/Ay) - tanh{{y - y2)/Ay)]/2 
with yi— km, y2= 60. kms, Ay =100 me- 
ters and a Webster profile along z given by 
An atan[12.63(z — ho)/ Ah]/!: with ho ~ 175 
meters. Ah —100 meters and An= -20.0 N units. 
The launching angles are 0.1, 0.2, 0.3, 0.4 and 
0.5 degrees in the vertical plane with and 0.001 
degrees in the horizontal plane. The TX is at 125 
meters along the z axis and the TX-RX antenna 
separation is 60 kms. 



